Lattice-gas model for collective biological motion 



Zoltan Csahok and Tamas Vicsek 

A simple self-driven lattice-gas model for collective biological motion is introduced. We find 
weakly first order phase transition from individual random walks to collective migration. A mean- 
field theory is presented to support the numerical results. 
PACS: 05.50.+q, 64.70.-p, 05.60,+w 

I. INTRODUCTION 

One of the most interesting aspects of evolution is the emergence of multicellular organisms due to the appearance 
of cooperation and differentiation of eucariotes. Although much research has been done along this line, some of the 
related basic questions are still open. As a natural step towards the understanding of the physical and physico- 
chemical background of self-organization of microorganisms several authors considered relatively simple systems such 
as the development of bacterial colonies. 

The growth of bacterial colonies having complex geometries has been extensively studied recently Jl]-^]. One of 
the mean aspects was the fractality |(| of the growing colonies. It has been found that the framework of diffusion 
limited growth fits well these phenomena In addition, the various morphologies of growing colonies have been 

experimentally investigated recently Q and a dynamical model has been introduced which incorporates a wide range 
of effects relevant to the phenomenon of collective bacterial growth and motion, for example chemotaxis. In a 
further related model Ym, aimed at describing the collective motion of self-driven particles (such as bacteria), a quasi- 
continuous approximation has been used with rules (particles moving with the same absolute velocity take on the 
average direction of motion of the neighboring particles) based on biological assumptions. As a main result it has 
been shown that spontaneous breaking of rotational symmetry can occur as the density of particles is increased or 
the level of random noise (i.e., the temperature) is lowered. The transition has been found to be continuous. 

One of the basic differences between living and azoic systems is that living organisms are self-driven: they can 
transform energy gained from food into mechanical energy which allows them to change their position. As the 
simplest example we can take bacteria [^|^], having various ways of motion. One of the mechanisms is motion by the 
means of flagella: the bacteria have flagella functionally analogous to a propeller attached to a motor. The motion 
of organisms is not under control of an external field, as is common in physical systems. Instead, the environmental 
effects cause only a change in the local velocity of the organisms. Since living objects are capable of communicating 
in various ways (ranging from the sensing of chemicals to verbal communication among humans), an organism is in 
continuous interaction not only with its environment but also with other organisms in its neighborhood. Thus, in the 
first approximation a system of organisms can be considered as an open interacting multi-particle physical system. 
Then, one can attempt to apply the methods recently developed in the investigations of complex systems [ [UMl| . 

In this paper we present a simple lattice gas model for the collective motion of self-driven particles. Similar approach 
has been applied to traffic systems [[l2] [jdj which also belong to the class of self-driven systems. Further approaches to 
self-driven systems include reaction-diffusion description [|15| , investigation of the related integrodifferential equations 
fusf , molecular dynamics jL7| and cellular automata [fl8| , [19| . 

The aim of this paper is to extend the usual statistical physical description for a particular case of collective motion 
in systems of living objects. First we introduce our model, then we give theoretical description of the problem. In 
Section IV we present the numerical results and in Section ^ we summarize our results. 



II. THE MODEL 

Our model is defined on triangular lattice of L 2 sites with unit lattice spacing and periodic boundary conditions. 
We put N particles (bacteria) on the lattice, where N is not necessarily smaller than the number of lattice sites. The 
density of the particles is defined as 

Each site can be either empty or occupied by one or more particles. If more than one particle is present at a site 
then in the calculations only the lowest one will be considered, where the lowest particle is defined as having the 
smallest random number previously assigned to each particle. 
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The particles are characterized by their position and velocity v, (i = 1 . . . N) which is of unit length (|vi = 1) 
and can point in any of the lattice directions (u a , Fig. 11.). 

At one time step positions and velocities of all particles are updated simultaneously according to the following rules: 

1. (a) for the particles which are not the lowest at their position we assign a random direction; 

(b) for the particles which are the lowest at their site we choose a new velocity u Q from a Boltzmann distribu- 
tion; 

P{u a ) = -|exp(-/3u Q ^2 Vj), 

where Z is a normalizing factor so that X)a=i-P( u a) = 1; an d P is ^/T (fcs = 1). The summation goes 
over the nearest neighbors which are in lowest position (Inn). 

2. every particle is moved one lattice unit in direction of its velocity: 



Note that the last step may result in sites with occupancy higher than one, this is the reason why we have to 
deal with such cases. The motivation for Step 1(a) is that we try to minimize the effect of multiple occupancy by 
letting the extra particles to diffuse away. The temperature parameter is not connected to the ambient temperature 
of the bacterial colony, it is rather an effective value which depends on many external parameters, as for example 
food concentration and agar humidity. The case of high food concentration is likely to be represented by high T 
values since then the bacteria can move faster and do not need coordinated behavior to extract food from the agar. 
On the other hand, when there is a food shortage the bacteria tend to cooperate, which results in a lower effective 
temperature. The above model is in its spirit close to the continuum model of self-driven particles [Q, however, the 
present version has a number of new features which had to be introduced because of its discrete nature. 

One of the quantities of interest is the average velocity of the particles which we shall consider as the order parameter 
and define as 
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N 



(2) 



Obviously < m < 1 holds. To have a closer analogy with spin systems we define a Hamiltonian 



H -- 

2 

in accord with the simulation rule Step 1. Those particles which are not the lowest at their position do not give 
contribution to the energy, they are regarded as a free gas. The energy per particle is 

/ \ E 
{£) = N- 

Having introduced the energy it is straightforward to define the heat capacity per particle 

de 
C= df- 

Although we have similarities with spin systems our model differs in a very specific way: the spins in our model are 
moving and this spatial dynamics is coupled to the spin dynamics. 

Fig. 11. shows a possible time evolution of the position and velocities of five particles. The particles are lettered 
from A to E and the arrows show the direction of their velocity (v^). At time step b) they form a cluster (containing 
a doubly occupied site) which then gradually breaks up. 
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III. THEORY 



Our system in closely related to the 6-state Potts model [^0| since we have q = 6 possible velocity states for a 
particle. Unlike the Potts model these states are not orthogonal, and we have an essentially non-equilibrium system, 
nevertheless, a mean field theory can be constructed in a similar way pl[ |. 

First we introduce a mean field Hamiltonian instead of Eq. (||) 

where the summation goes over all lowest (1) particles, not only the nearest neighbors. The mean field energy function 
can be written as 

1 6 

EmF = --^Ng cS 2J XaUayX^, (5) 
a, 7—1 

where g c g — 1 — exp(— g) is the effective density, i.e., a site has on average 6g e s occupied neighboring sites, x a is the 
fraction of particles travelling in the lattice direction a QZ Q =i x a — 1) and \J ai — u a u 7 which for the case of the 
Potts model would be simply U ai = 5 ai . 
The average energy per particle is 



Emf 



= ~\q*& X/ XaUajXj. (6) 

The entropy per particle is 



N 2 L 

a, 7=1 



smf = — y^^a In ate, 



a=l 

so for the free energy per particle one gets 



6 / 1 6 \ 

^ ^ I x a In x a ^^eff/^^ ^ ] U a ^Xj j 
Q = l \ ~ 7 = 1 / 



We intend to find the configuration x a which minimizes the free energy /mf- Since all the lattice directions are 
equivalent we can look for a solution in the form of 

1 5 

x\ = - + -m M F (8) 
6 o 



and 



1 m MF , n , 

x a >i = — , (9) 

6 6 



where mMF is the mean field order parameter which satisfies the relation 

6 



a = l 



according to Eq. (|2|). Substituting Eq. (||) and Eq. (||) to Eq. (^) after a bit of algebra one gets for the mean field free 
energy 

a ± 1 a 2 i 1 , 5(1-to M f) 1 1-tomf 
P/mf = -^g c «P m MF - log - + log 

1 + 5TOMF, l + 5m M F , lrA 

+ a l0 S a ■ ( 10 ) 
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For high temperatures (T > T c ) this function has its minimum at uimf = which corresponds to the disordered 
state of the system. At the critical temperature, which can be derived from /mf and in our case it is 

c 3.353' [ ' 

a non-trivial minimum appears. The phase transition, like in the 6-state Potts model |22|| , is first order. The jump in 
the order parameter in this approximation is given exactly by 

AmMF = 0.8. 



IV. NUMERICAL RESULTS 



We have studied our systems by Monte-Carlo simulations. For initial configuration we chose random distribution 
for the position and velocity of the particles. Typical configurations of the system for various temperatures T and 
particle densities g are shown in Fig. ??. It can be easily seen that at low temperature the particles tend to form 
clusters as it can be expected. 

We have performed several long-time runs to obtain the behavior of the quantities defined in Section [□]. as a 
function of T and g. We used various system sizes (L) up to 40 for high densities and up to 200 for low densities. The 
limiting factor was the convergence time which for the case of our largest system was in order of 10 6 sweeps of the 
system. Fig. ??. demonstrates the order parameter as a function of the temperature for g = 0.9. The estimated jump 
at T c is smaller but close to the value obtained by the mean field theory. In Fig. ??. we present the average energy 
which is also subject to a finite jump at T c . These two figures suggest that a first order phase transition takes place 
at T — T c in agreement with the theoretical prediction. A strong evidence supporting this idea is presented in Fig. 
??. where we have plotted the distribution P(e) of the energy values for a number of different temperatures below 
and above T c . One can clearly see a gap in the distributions at intermediate energies which is a unique feature of first 
order phase transitions Q . The inset in the figure shows the distribution for Ising type interaction of non-moving 
particles in the same system where the transition is known to be second order. In Fig. ??. we present the temperature 
dependence of the heat capacity which is the measure of the broadness of the energy distribution. A characteristic 
peak can be observed near T c . The position of the peak is shifted for various lattice sizes due to finite size effects. 

We studied the behavior of the model also as a function of density of particles (g). Fig. ??. shows the temperature 
dependence of the average energy for various densities obtained. The transition is present even for very small densities 
although the position of the critical temperature lowers. In Fig. ??. we have plotted the dependence of the transition 
temperature on the density. There is a natural distinction between the high and low density regimes of the system: 
at the percolation threshold the behavior of the system is expected to change. In fact we observe a change in the 
dependency of the critical temperature below the percolation threshold of the triangular lattice (g = 1/2 and g e g w 
0.39) at g c ff « 0.25 which corresponds to density g w 0.29. The values of the measured critical temperatures are higher 
than the one obtained from Eq. ( |Tl| ) which shows the boundaries of applicability of our mean-field approximation. 

The behavior of the average energy of the Potts model near its transition temperature can be characterized by two 
exponents a^~^ and aW [ p2[ . These exponents are present due to the weakly first order nature of the transition. The 
temperature dependence of the average energy is given by 

(e) = £ (-)-A(-)(l-T/T c ) 1 - Q< " ) (12) 

for T < T c and similarly 

(e) =e {+) +A {+ \l-T c /T) 1 - aW (13) 
for T > T c . The difference between and e^ - -* is equal to the energy jump during the phase transition. In Fig. ?? 



and Fig. ?? we plotted the energy differences versus the temperature according to Eq. (12) and Eq. ( |l3| ) for g = 0.9. 
The exponents obtained from the slopes are 

a (+) w 1 - 0.73 = 0.27 

and 

a { - ] w 1-0.5 = 0.5. 

These values are different both from those of the q = 6 state Potts model (a^ w 0.7 and er ^ « 0.7) and in part 
different from the corresponding mean field values (a^ = = 1/2). 
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V. CONCLUSION 



We have presented a lattice-gas model for collective biological motion. We have shown both numerically and 
theoretically that weakly first order phase transition takes place in our system separating the phase with zero net 
transport and the ordered phase with non-zero average velocity. We find the exponents and a^ - -* are differ from 
the ones of the q = 6 Potts model and from mean field values. This difference can be attributed to the fact that 
although we have similarities with spin systems our model differs in a very specific way: the spins in our model are 
moving and this spatial dynamics is coupled to the spin dynamics. It is remarkable that the behavior of the present 
lattice model is qualitatively different from that of the analogous continuum model M. While in the continuum 
model and in a directly related continuum equation for a two-dimensional dynamic XY model p4[ a second order 
transition was observed, in our lattice gas version the transition is more complex and has a first order component. 
Such discrepancies, however, are not unfamiliar even in two-dimensional equilibrium systems: in particular, there is 
no long range ordering in the equilibrium XY model p5|| having continuous symmetry, while its discrete counterparts 
(i.e., the Ising model) exhibit second order phase transition. 
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VII. FIGURES 



Fig.l. One lattice site, the six lattice directions are shown. 

Fig. 2. (a-d) Possible time evolution of our model (a — > d) shown on a small portion of the lattice. Note the double 
occupancy in time step b. 

Fig. 3. (a-d) Some typical snapshots of the system at (a) high temperature, (b) intermediate temperature, (c) low 
temperature at g = 0.9 and (c) intermediate temperature at a lower density. Note the appearance of ordered clusters. 
(Only the particles in lowest position are drawn.) 

Fig. 4. The order parameter (m) as a function of the temperature (T) for density g — 0.9. 

Fig. 5. The average energy ((e)) versus the temperature in the same systems as on the previous figure. 

Fig. 6. The energy distribution (P(s)) for various temperatures below and above the transition (g = 0.9). Note the 
energy gap between e w — 1 and e « —3.5. Inset shows the distribution for Ising spins instead of mobile particles 
where the transition is continuous. 

Fig. 7. The heat capacity (c) versus temperature graph for the systems as on Fig. 4. The dotted lines are guide to 
the eye. 

Fig. 8. The average energy as a function of temperature for various densities of particles. 
Fig. 9. The critical temperature as a function of the density of particles. 
Fig. 10. Energy difference versus temperature for g = 0.9. (T c = 0.413, = —0.76) 
Fig. 11. Energy difference versus temperature for g = 0.9. (T c = 0.413, = —3.75) 
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FIG. 1. 



FIG. 2. 
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FIG. 3(b). 
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